Fourier transformation and response functions 
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We improve on Fourier transforms (FT) between imaginary time r and imaginary frequency cj n 
used in certain quantum cluster approaches using the Hirsch-Fye method. The asymptotic behavior 
of the electron Green's function can be improved by using a "sumrule" boundary condition for a 
spline. For response functions a two-dimensional FT of a singular function is required. We show 
how this can be done efficiently by splitting off a one-dimensional part containing the singularity 
and by performing a semi-analytical FT for the remaining more innocent two-dimensional part. 



Quantum cluster cluster theories, such as the dynam- 
ical cluster approximation (DCA) or the cellular dy- 
namical mean-field theory (CDMPT)ji make it possible 
to calculate dynamical quantities, e.g., electron Green's 
functions or response functions, for strongly correlated 
systems. These calculations, however, are numerically 
very demanding. In the Hirsch-Fye 2 method for solv- 
ing the resulting cluster problem, one has to switch be- 
tween imaginary times r and imaginary frequencies to n , 
which requires Fourier transforms (FT). In this paper we 
address efficient methods for performing FT in this con- 
text. In the weak-coupling version 3 -!^ of continuous time 
approaches 3-6 properties can be very easily measured di- 
rectly in frequency space and no FT is needed. 

The FT of the electron Green's function has to be per- 
formed with great care, and it is important to obtain the 
correct asymptotic behavior for large u) n £~— Inaccura- 
cies for large u> n can lead to problems, for instance in the 
DMFT description of a Mott transition.- The large ui n 
behavior can be greatly improved if exact moments are 
calculated. The FT can then be performed using a nat- 
ural spline, which works well for a symmetric half-filled 
Hubbard model.- For a nonsymmetric or doped Hubbard 
model, however, the accuracy is not optimum. Here we 
show how this can be improved by introducing "sumrule" 
boundary conditions for the spline interpolation. 

The main part of this paper deals with response func- 
tions. This is important not only for studying suscepti- 
bilities but also for diagrammatic extensions of the dy- 
namical mean-field theory, such as the dynamical ver- 
tex approximation^ and the dual fermion method.— The 
most convenient way of of calculating them within DCA 
is to work directly in w n -space. Yet, in the Hirsch-Fye 
method this requires FT for functions gfc, Tj) depending 
on two imaginary times, Ti and Tj, for each Monte- Carlo 
step. The FT is complicated by the fact that these func- 
tions have singularities and that the precise asymptotic 
behavior is difficult to work out. We show that a very 
efficient solution is to split up g{ji,Tj) into two parts. 
One part, g^iji — Tj), depends only on the difference 
Ti — Tj and contains the singularity. The second part, 
Sg(Ti, Tj) depends on both variables independently but is 
well behaved. To perform spline interpolations in both 
variables is very time consuming. We show how the FT 
can be performed in a very efficient way for Sg(ri, Tj). 



If G(t) is only known for discrete values of r sepa- 
rated by At, a direct FT can only give accurate results 
up to LJ n ~ 1/At. To improve the accuracy for large 
ui n , on can subtract a model Green's function, G m (r), 
from G(t), where G m has the right asymptotic behavior. 
The difference, G(t) — G m (r), is then FT and G m (iuj n ) 
is added. For instance, G m (r) can be obtained from per- 
turbation theory 13 Alternatively, one can calculate the 
lowest moments of the spectral function exactly from ap- 
propriate expectation values^ G m (iuj n ) can be chosen so 
that these moments are exactly satisfied, meaning that 
the corresponding coefficients in a (l/w„) expansion of 
G m (iui n ) are correct. G(r) — G m (r) is then FT using 
a natural spline^ If the FT of G(r) — G m (r) correctly 
gives the lowest moments equal to zero, it follows that 
the corresponding moments of G(iuj n ) are correct. 

To illustrate this, we consider the relation between the 
spectral function, A(uS), and G(r) 



G(t) = 



1 



(i) 



where f3 = 1/T and T is the temperature. This gives the 
following, very useful, relations 



G (n) (0) + G (n) (J3) 



{-u) n A{u)duj = (~l) n M n 



(2) 

where G^^t) = d"G(r)/dr n and M n is the nth mo- 
ment. For large oj n , G(u n ) ~ J2k Mk/ {io*> n ) k+1 ■ For the 
difference AG(r) = G(t) — G m (r), it then follows that 

AG (n) (0) + AG (n) 03) = n = 0,l,2, (3) 

if G m has the correct Oth, 1st and 2nd moments. 

Let AG(t) be given for n r-values, t\, r„. In a 
third order spline, each interval between two r-values is 
interpolated with a third order polynomial. It is required 
that the polynomials give the correct AG(Ti) and that 
the first two derivatives are continuous at each r^. There 
are then 4(n — 1) unknown coefficients and 4(rt — 2) + 2 
conditions. We then have to provide two more conditions. 

For a natural spline, it is assumed that AG^ 2 ^(0) = 
AG( 2 >(f3) = 0, while the first derivatives are left open. 
The assumption about AG 1 - 2 ' is too strong, since we 
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FIG. 1: (color on-line) Deviation of the local £(2) (z = iuj n ) 
from the asymptotic form E(z) ~ a + b/z for natural and 
sumrule boundary conditions. For natural boundary condi- 
tions the frequency independent part a has an error. 

only know that AG^ (0) + AG (2 ^ {(3) = 0. The sec- 
ond moment is nevertheless correct, since it only de- 
pends on the sum. The incorrect assumption about 
AG^- 2 \ however, influences the estimates of AG' 1 -* and 
AG (1) (0) + AG (1) (/3) = is in general not satisfied. This 
leads to incorrect results already for the important first 
moment. Due to symmetry, the natural spline may give 
AGW(O) + AGW(/3) = in special cases, e.g., at half- 
filling for the symmetric Hubbard model. 

A better approach is to use the two conditions of 
Eq. @ for n = 1 and 2. We refer to this as the sumrule 
spline. This gives correct first and second moments and 
large uj n behavior, even if the estimates of AG^ (0) and 
AGW(/3) individually are not accurate. 

FigfT] compares the two approaches for the two- 
dimensional (2d) Hubbard model in the DCA. We have 
used t — —0.4 and t p = — 0.3i for the nearest and next 
nearest neighbor hopping, respectively, U — 8|i| for the 
on-site Coulomb interaction and f3\t\ = 12. All energies 
are in eV. The occupancy is n = 0.9 and we consid- 
ered an eight site cluster. The number of r-points is 
N T = 120. For large ui n , the local Green's function be- 
haves as G(z) ~ l/(z — a — b/z), where z = iuj n and 
a = Mi and b = M 2 — M 2 are given by the first two 
moments. If we define S(z) = z — G _1 (z), we have 

Y,(z)~a + b/z (4) 

for large \z\. The figure shows that Re S(z) — } a is 
obtained for the sumrule spline but not for the natural 
spline. Also for Im £(z) ~ b/z (not shown in the fig- 
ure) the natural spline is less accurate than the sumrule 
spline, since an error in Mi also enters in b, but for large 
I z I the error for Im E(z) is much smaller than for Re £(z). 
As the write up of this work was being finished, we be- 
came aware of very similar approaches in the thesises by 
Comanao^ and by GullJ^ 

We now discuss the FT for response functions, and use 
DCA as an illustration. The response function is calcu- 
lated for a cluster in a bath. From the cluster response 



function vertex corrections are deduced. The Bethe- 
Salpeter equation for the lattice is then solved, assuming 
that the vertex corrections are the same as for the clus- 
ter. The embedded cluster problem is solved for imag- 
inary time, while the Bethe-Salpeter equation is solved 
for imaginary frequency. The necessary FT to imaginary 
frequencies is numerically difficult. 

We consider the electron-hole response function 

1 [P ^ [P [P 
X <T<T '(q,k,k) = --^ / dn / dr 2 / dr 3 / dr 4 

P Jo Jo Jo Jo 
xexp{-i[uj n Ti - (w n + v)t 2 + (uy + v)t 3 - u; rl /T 4 ]} 

x (Tr [<L (ri )c k+q . (r 2 )c[ , (r 3 )c k , (r 4 )]> , (5) 

with the compact notations k — (k, lo„) and q — (q, v). 
Here T T is a time ordering symbol, creates an 

electron with wave vector k and spin a and c(r) = 
exp(i5/r)cexp(— Ht), where H is the Hamiltonian. 

The c operators are contracted pairwise and their ex- 
pectation values are calculated. This is done for a very 
large number of configurations. It is convenient to per- 
form the FT to imaginary frequency and reciprocal space 
for each configuration and to store the results in the k- 
and q- variables rather than storing the results in imag- 
inary time and real space. We then need 

g^fcVEE E e - iKn ~^' Tjl (6) 

i=l j=l RiRj 

xe i t k - R - k '- R ^(T T [ C%CT (r J -)4 iCT (r i )])(Ar) 2 , 

where the integrals over r have been replaced by sums 
over N T discrete values of r separated by At and R is 
a site index. Although this has the form of a Green's 
function, it is calculated for a specific configuration and 
only the zeroth moment is known. This makes it harder 
to perform a FT. There is a singularity at Ti — Tj, which 
makes a straightforward spline in and Tj less useful. 
The singularity can be handled by treating 7j < Tj and 
Tj > Tj separately. But a spline in two variables is still 
numerically very demanding because of the large number 
of points needed. To see this, we simplify the calculation 
in Eq. ([6]), by splitting it in two parts. Thus we calculate 

f a (k,R J ,T J ) = J2e- l ^ T '- k - R ' ) (7) 

and 

g a (k,k) = E e i K"-i- k '-^)/ (r (fe ) R J ,r J -). (8) 
This gives 

Xaa'^k^k) = -jpl9a(k + q,k)g a >(k' ,k + q) 
-g„(k' ,k)g„{k + q,k' + q)5 aa ,}(AT) 2 (9) 
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Let the number of sites be N c , the number of w n - values 
N u and the number of ^-values N v . The number of k- 
values is then also N c . Let N% be the number of r-points 
for which the correlation function in Eq. (J7J is known 
and N* the number of r- values needed to obtain an ac- 
curate FT. Eq. (0 and Eq. ((8]) then require of the order 
of 2{N U + N„)NZN°N* and 2(jV w + N„) 2 N°N* opera- 
tions, respectively, for each configuration. Here we have 
assumed that the spline in the first t- variable is only done 
for each of the N£ values of the second r-variable. After 
the corresponding FT has been performed, the second 
variable is splined and FT. The calculations can easily 
be arranged so that the time needed for calculating the 
exponents is negligible and efficient machine routines can 
be used for the multiplications. Still, the calculations are 
very time consuming if is large enough to give accu- 
rate FT. We therefore follow a different route, reducing 
the time requirement for Eqs ([JJ |5]) very substantially 
and requiring no interpolation of the r-variables. 
We first notice that 

ffo-Ri ,R 3 - (ji > Tj ) = ([T T C R . a (T j )4 li(7 (Ti))]) (10) 

depends on Tj and Tj individually and not only on their 
difference, since it is calculated for one particular config- 
uration. However, we can separate it as 

5<tr, Rj(n, ^ ) = 5° Ri Rj (n - ^ ) + SgrB* Rj (n,Tj), ( 1 1 ) 

where 



3°RiRj ( T i) 



1 

iVV 



N T 



g^njRf {n+j-i , Tj 

i=i 



(12) 



only depends on one r-variable and we have defined 

9*-Ritlj(n+j-l,Tj) = -gaH t -R. 3 {Ti+j-l-N T ,Tj) if i + j - 
1 > N T Or gvUtRj (n+j-l,Tj) = — gaUiRjin+j-l+NriTj) 

i + j — 1 < 1. Here, fl^R r ( t * s n0 ^ a noninteracting 
Green's function, but the time translationally invariant 
part of go-RiR,- (t,-_(-j_i, Tj). The singularities are now in 
5° R . Rj (r), and Sg^R^ (n, Tj) is free of singularities, and 
can more easily be Fourier transformed. 

Since 5° R . Rj . (r) only depends on one variable, it can 
easily be FT using a spline. Alternatively, we can use 
Filon's rule^ where second order polynomials are fit- 
ted to the -/V T r-points. These polynomials are then FT 
analytically. Even for w„At 3> 1, the FT can be very ac- 
curate. This automatically gives the appropriate \/(iu) n ) 
behavior for large ui n , due to end point corrections. 

It is possible to FT Sg a -R i R.(ri,Tj) by performing a 
Filon's rule for first Tj and then for Tj. However, we have 
found it preferable to fit a two-dimensional polynomial 



aoo + o-iqt + <iqit' + 



o,\\tt 



(13) 



to the values of Sg^R^ (t, t') in the points (t^t, ), 
(t 4+ i,t,-), (n,Tj +1 ) and (n + i,Tj + i). This is multi- 
plied by the appropriate exponent and integrated an- 
alytically. Substantial simplification follow from the 




FIG. 2: (color on-line) The relative accuracy of the FT of 
Eq. (|16[) according to the approach of Eq. (|14|l (Filon) or 
using the trapezoidal rule (Trapez). 



fact that <5<7 CTRiRi (t, t ) is antiperiodic and exp[z(w n T 
w n' r ']^fffRiRj ( r ' T ) ^ s P er i°dic in r and t'. Then 



rP 
dr I 
Jo 



JcK,^,At/2) £ ^ 7 



r,r j (t,t)= (14) 
»' T,, ^R,R i (r„r J ), 



where n = (i - 1)At and (N T + 1)At = /3. Here 
c(x,y,A) = 

e -iA(x+y)[ b x b -y 

e iA ^ +y) [b x bo y + 6§&i 



h x h~ y 
°0°l 



&5V 



where 



°o 



1 1 

(e 

x 



ixA 



h x h~ y 
°0°1 



-ixA\ 



&5V 

ft?6 -» 



(15) 



Ax A 



(1 - ia;A) 



-ixA 



(1 + isA)] 



This approach can easily be extended to the case of a 
nonuniform grid. 

If ^ffo-RiR,- (t, t ) were a very smooth function, a more 
accurate integration method could be devised by fitting a 
polynomial of higher order. However, since 8g is obtained 
for a specific configuration, this does not seem useful. 

To test the method, we have FT a function 



f(n,T 2 ) =^Oij( 



p/2 



/3/2 )i( T2 



0/2, 



0/2 



(16) 



where is only nonzero for odd values of i and j to as- 
sure that the function is antiperiodic. Specifically, we 
chose an = 0.7, 013 = 1.3, ais = 0.9, 031 = —1.2, 
033 = 1-5, a 35 = -0.6, a 5 i = -0.8, a 53 = 1.1, a 55 

0.7. We used the frequencies u> n — 13.5(27r//3) and 
LJ n = 9.5(27r//3), where /3 = 15. Fig. [5] shows results ob- 
tained by using Eq. (|14l) (Filon) or the simple trapezoidal 
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FIG. 3: (color on-line)The quantity i/ 2 II(j/) as a function of 
v for different values of N T . The figure also shows results 
when g° has been split off [Eq. (|f ![) ] but the trapezoidal rule 
was used for 8g (g°) or the total g was integrated using a 
trapezoidal rule (Trapez). 

rule (Trapez). In the figure, the approach of Filon leads 
to a comparable accuracy as the trapezoidal rule for a 
N T that is almost one order of magnitude smaller. 

In this Filon like approach the exponent is treated ex- 
actly and the error in the FT is entirely due to the limited 
information about the function to be FT. It is then no 
gain in adding extra points by interpolating the function 
to be FT. In a Hirsch-Fye approach this means that we 
put N T — the number of points determined by the 
discretization used. 

From Xaa'iiQ'k'k ) [Eq- ©] we can calculate II(r) = 
(j(r) • j(0))/(3iV), where N is the number of lattice sites 
and j is the current operator. The FT of H{t) is related 
to the optical conductivity er(w) via 

1 foo 2 

n M = ~/_^.(.)dc, (it) 

Eq. (JTTJ) shows that v 2 H(v) approaches a constant for 
large v. Problems of the FT should show up in particular 
for large v and the accuracy should increase with N T . 



We then choose N T so large that v 2 Ii(v) is constant for 
large values of i/- values. This should then be an accurate 
result. 

Fig. [3] shows results for v 2 W{v) for the 2d Hubbard 
model. The parameters are the same as in Fig. [TJ except 
that f3 = 15. The bath obtained for N T = 160 was used 
also for N T — 60. For N T = 160, v 2 Tl{v) is constant for 
large v over the whole range shown. The comparison with 
N T = 60 suggests that the FT is quite accurate at least 
for vAt < 2 and it stays fairly accurate for substantially 
larger values i/Ar. The deviation between N T = 60 and 
160 could also be due to other inaccuracies for N T = 60 
than the FT, and in that case the FT is accurate for even 
larger is At. The figure also shows a calculation where we 
split off g°, [Eq. dTTJl], and FT it using Filon's rule, but 
FT 5g using the trapezoidal rule (g° in the figure). We 
also performed the FT on the full <?, without splitting 
off g , using the trapezoidal rule (Trapez in the figure). 
The figure shows that for iV£ = 60 both approaches fail 
dramatically for large v. 

To summarize, the FT of the Green's function can be 
improved by using a spline with sumrule boundary con- 
ditions. This gives a G(iu> n ) with correct first and second 
moments, while a natural spline in general gives an in- 
correct first moment. To calculate a response function, 
we need a FT a function g(ji,Tj) with a singularity. We 
show how a g° can be split off, which only depends on the 
difference Tj — Tj and which contains the singularity. This 
function can be FT very accurately. For the rest, Sg, we 
developed a two-dimensional FT in the spirit of Filon's 
rule. This leads to accurate results, even if g{ri,Tj) is 
only known on a rather sparse mesh. "After this paper 
had been submitted, an alternative prescription for effi- 
cient Fourier transforms of two-particle Green's functions 
has been proposed by Kunes^ 
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